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Resolution  Requirements  for  the  Numerical  Computation  of  Tonal 
Noise  in  Compressors  and  Turbines  of  Aeroengines 


T.  Hiittl1,  G.  Kahl,  F.  Kennepohl  and  K.  Heinig 

MTU  Aero  Engines  GmbH 
Dachauer  Strasse  665 
D-80995  Munchen,  Germany 


Abstract 

The  time  linearized  Euler  method  Lin3D  is  applied  to  two  sets  of  test  cases.  2D  wave  propagation  test  cases  in 
homogeneous  flow  are  used  to  quantify  numerical  dissipation  and  dispersion  of  the  discretization  scheme.  The 
minimum  number  of  mesh  diagonals  between  two  wave  fronts  has  been  found  to  be  an  appropriate  measure  of  the 
resolution  of  a wave.  Correlations  have  been  found  that  characterize  the  dissipation  and  dispersion  behavior  of  the 
code  and  therefore  the  resolution  requirements  for  a given  flow  simulation.  The  transmission  and  reflection  of  plane 
sound  waves  incident  upon  a single  cascade  of  finite  plates  has  also  been  calculated  with  Lin3D  and  compared  with 
an  analytical  method  of  Koch  [5],  The  computed  ratios  of  transmitted  or  reflected  to  incident  pressure  wave 
amplitude  agree  well  with  the  analytical  solution,  even  for  scattered  modes. 


1.  Introduction 

The  noise  produced  by  jet  engines  has  become  one  of  the  major  concerns  in  today's  pollution- 
and  noise-conscious  society.  The  jet-exhaust  noise  has  been  reduced  significantly  with  the 
increasing  by-pass  ratios  of  modem  turbofan  engines.  But  now,  the  pure  tonal  noise,  produced 
by  rotor-stator  interaction  in  compressors  and  turbines  can  become  a dominant  source  of  engine 
noise.  Together  with  the  noise  generation  and  radiation,  the  sound  transmission  through  blade 
rows  in  axial-flow  turbomachines  is  also  of  great  importance.  Various  calculation  methods 
ranging  from  analytic  methods  for  simplified  flow  and  geometry,  over  empirical  methods  based 
on  correlations  calibrated  by  experimental  data,  to  the  wide  range  of  numerical  methods  have 
been  established  within  the  aeroacoustic  design  process  of  turbomachines. 

In  this  paper,  we  focus  on  a numerical  calculation  method,  based  on  the  time  linearized  Euler 
equations.  This  CFD  code  Lin3D  developed  at  MTU  is  already  used  for  aeroelastic  design  of 
aeroengine  compressors  and  turbines.  It  may  also  be  applied  to  the  aeroacoustic  design  of  these 
engine  components.  In  contrast  with  analytical  methods  a more  realistic  description  of  blade  and 
duct  geometry  is  possible.  Compared  with  Navier-Stokes  methods,  the  computational 
requirements  of  Lin3D  are  much  smaller.  Before  applying  a code  to  a new  kind  of  problems,  its 
requirements  and  weaknesses  have  to  be  studied.  It  is  the  aim  of  the  present  study  to  investigate 
the  resolution  requirements  of  Lin3D  for  aeroacoustic  computations. 

First,  the  wave  propagation  in  homogeneous  flow  has  been  computed  in  order  to  quantify  the 
effect  of  numerical  dissipation  and  dispersion  on  the  solution.  Then,  the  problem  of  the 
transmission  of  a sound  wave  through  a blade  row  has  been  solved  where  the  accuracy  of  the 
code  can  be  benchmarked  with  an  analytical  method. 
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2.  Numerical  Method 

The  CFD  code  Lin3D  is  a time-linearized  Euler  method  for  3D  cascade  flows  [2,3,4].  The 
unsteady  flow  is  assumed  to  be  a small  harmonic  perturbation  of  the  non-linear  steady  flow,  so 
that  the  steady  flow  problem  is  decoupled  from  the  unsteady  problem.  As  long  as  the  wave 
amplitudes  remain  moderate,  the  higher  order  terms  in  the  governing  equations  derived  under 
this  assumption  can  be  neglected  and  the  describing  unsteady  flow  equations  become  linear.  The 
small-disturbance  Euler  equations  are  transformed  into  the  frequency  domain,  so  that  the  results 
are  given  as  amplitudes  and  phases  of  the  flow  variables  in  the  whole  (3D)  computational 
domain. 

The  base  flow  for  the  time  linearized  Euler  code  has  been  prescribed  analytically  for  the 
investigations  discussed  here,  but  it  can  be  the  result  of  a steady  3D  Euler  computation,  as  well. 
Lin3D  is  restricted  to  small  disturbances  of  the  mean,  inviscid,  steady  flow  and  its  applicability 
is  limited  to  cases  with  (predominantly)  attached  flow.  The  code  is  based  on  the  linearization  of 
the  steady  Euler  method  described  in  [1]  and  is  an  extension  of  the  code  described  in  [2,4], 

Lin3D  uses  a Finite-Volume  scheme  with  a 3-step  Runge-Kutta  time  integration  method.  A 
structured  H-type  mesh  with  node  centered  variable  arrangement  is  required  to  describe  the 
geometry  of  the  flow  domain. 

Lin3D  has  proven  to  be  sufficient  for  most  cases  concerned  with  turbines.  It  is  a fast  and  robust 
tool  routinely  used  in  the  design  process  to  assess  flutter  stability  and  forced  response  of 
compressor  and  turbine  bladings.  It  features  deforming  meshes,  non-reflecting  boundary 
conditions,  arbitrary  eigenmodes  (including  chordwise  bending),  forced-response  calculations 
for  generalized  forces  due  to  up-  and  downstream  disturbances  and  the  possibility  of  calculation 
of  acoustic  modes.  Lin3D  runs  on  a wide  variety  of  platforms  ranging  from  workstations  to 
multiprocessor  vector  supercomputers,  with  the  option  of  parallel  computing  to  speed  up  the 
turnaround  times. 


3.  Wave  propagation  in  homogeneous  flow 

3.1  Computational  domain  and  setup 

The  accuracy  of  the  code  for  wave  propagation  in  homogeneous  2D  flow  has  been  tested  for 
numerous  upstream  and  downstream  propagating  waves.  For  each  single  test  case  the  same 
equidistant  mesh  has  been  used,  where  the  vector  of  the  homogeneous  mean  flow  is  parallel  to 
the  mesh  lines,  see  Figure  3.1.  Physical  and  geometrical  parameters  of  the  cases  are  listed  in 
Table  3.1.  In  order  to  realize  a plane  2D  flow  with  the  3D  code  in  cylindrical  coordinates,  a 
small  ratio  of  channel  height  ( rlip  - rhuh  - 0.2 m ) to  tip  radius  ( rtjp  = 300 m ) has  been  chosen  and 

discretized  by  4 mesh  cells. 


Table  3.1:  Physical  and  geometrical  parameters  of  2D  wave  propagation  calculations. 
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Figure  3.1:  Flow  domain  of  2D  wave  propagation  calculations  (without  blades). 


Two  parameters  are  varied,  the  circular  frequency  co  and  the  incidence  angle  </>  of  the  pressure 
wave: 


_ _ _ /V*-* 

to A = 0-25—,  c»b=Y’ 

<p  = 0°,15o,30o,45°, — ,345° 


=*Y’  0)0  = 16_7r 

(not:  90°  ,270°). 


For  upstream  traveling  waves  (90° <0<27O°)  the  boundary  conditions  for  a pressure 
perturbation  have  to  be  specified  at  the  outlet  plane;  for  downstream  traveling  waves 
( 0°  < ^ < 90°  and  270°  < 0 < 360° ) at  the  inlet  plane.  All  waves  are  sinusoidal  in  the  y - 
direction.  The  amplitudes  are  sufficiently  small,  so  that  non-linear  effects  can  be  neglected  (peak 
to  mean  pressure  variation  of  0.1%  of  the  static  pressure). 


3.2  Resolution  of  the  waves 


Figure  3.3:  H-type  mesh  geometry. 
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Although  mean  flow,  flow  domain,  number  of  mesh  points  Nx , Nv  and  wave  number 


a 

are  chosen  to  be  identical  for  constant  co,  the  resolution  of  the  waves  is  not  equal  for  the 
simulation,  as  it  depends  on  the  incident  angle  <p . However,  the  wave  lengths  X also  depend  on 
<p  for  flow  at  Ma  * 0 . A characteristic  quantity  for  the  resolution  R in  discrete  domains  is  the 
minimum  number  of  mesh  cell  diagonals  between  two  wavefronts: 


R = min 


4 


X(a, ) is  the  distance  between  two  wavefronts  along  a line  under  an  angle  of  ai  to  the 
horizontal  axis,  see  Figure  3.2. 

X(a, ) = — . 

cos  or,  + tan  0 sin  or, 


Xx  can  be  obtained  from 

Xx  = , using  cx  = nMacos/?,  H — - — f-nMasin^  tan^  . 

CO  COS0 

Two  diagonals  and  d2  exist  in  a structured  H-type  mesh  (Figure  3.3).  Their  lengths  and  their 
angles  a, , a2  are  defined  by 

dx  = J A x2  + (Ax  tan  J3S  + Ay)2  , tan  a,  = tan  /3S  + — , 

Ax 

d2  = Jax2  + (Av  tan  - A yf  , tana2  = tan/?s , with  Av  = — , Ay . 

Ax  N x Nv 

3.3  Results 


Figure  3.4:  Unsteady  pressure  for  four  wave  propagation  test  cases  with  the  same  incident 
angle  <p  = 105°  and  different  circular  frequency  co : a)  (0A , b)  coB , c)  coc  , d)  coD . 
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Figure  3.5:  Unsteady  pressure  for  four  wave  propagation  test  cases  with  the  same  circular 
frequency  (Oc  and  different  incident  angle  0:  a)  0 = 165°,  b)  0 = 225°,  c) 
0 = 285° , d)  0 = 330°. 


A first  impression  of  the  wave  propagation  in  a 2D  flow  field  is  given  in  Figures  3.4  and  3.5. 
Figure  3.4  demonstrates  that  the  damping  of  the  waves  increases  for  higher  co  (smaller  k , lower 
R ).  In  Table  3.2,  the  values  of  the  resolution  R for  four  cases  with  the  same  circular  frequency 
(t)c,  plotted  in  Figure  3.5,  are  listed.  This  Figure  shows  that  the  damping  depends  on  the 
resolution  R , which  is  a function  of  the  incident  angle  0 of  the  wave. 


The  purpose  of  the  wave  propagation  test  cases  is  to  analyze,  how  strong  numerical  dissipation 
and  numerical  dispersion  can  affect  the  flow  solution.  Numerical  dissipation  and  dispersion  are 
unwelcome  side  effects  of  a discrete  solution  method  and  depend  on  the  discretization  scheme 
of  the  code.  Numerical  dissipation  leads  to  the  exponential  decrease  of  the  wave  amplitudes 
A(x)  in  (horizontal)  propagation  direction  and  can  be  quantified  by  the  decay  rate  Dec 

r A{x  + AxY' 


Dec  = —20  • log 


Dec  = —20  • log 


for  a downstream  propagating  wave 


for  an  upstream  propagating  wave 


Numerical  dispersion  leads  to  an  incorrect  propagation  velocity  of  a wave  or  to  a change  in  its 
propagation  direction.  Dispersion  can  be  quantified  by  the  phase  change  error  Pee,  the 
difference  between  analytical  Pa(x ) and  numerical  phase  Pn{x)  change 

Pee  = \{P„  (x  + Ax)  - Pn  (x))  -{Pa{x  + Ax)-Pa  (x))| 

Table  3.2:  Resolution  of  the  pressure  waves  of  the  four  test  cases  in  Figure  3.5. 


Case 

a 

b 

c 

d 

Incident  angle 

0 = 165° 

0 = 225° 

0 = 285° 

0 = 330° 

Resolution 

R=  22.4 

R = 3.92 

R = 5.43 

R = 40.0 
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Figure  3.6:  Decay  rate  Dec  (left)  and  phase  change  error  Pee  (right)  of  the  wave  propagation 
test  cases  drawn  over  the  resolution  R : A coA , O coB , coc , O coD , lines: 

Declq  = 9.2  • R~2 5 or  Pcelq  =1.1-  R~2 ! . 


Figure  3.7:  Decay  rate  Dec  (left)  and  phase  change  error  Pee  (right)  of  the  wave  propagation 
. test  cases  drawn  over  the  resolution  R for  two  different  convergence  criteria:  A 
low  convergence  criterium  ( Resp  =5-1 0-7 ),  O high  convergence  criterium 

( Res^  =5-1 0-5 ),  lines:  Declq  = 9.2  -R~25  or  Pcelq  = 1.1/T21. 

In  Figures  3.6  and  3.7,  the  decay  rates  and  the  phase  change  errors  of  the  wave  propagation  test 
cases  are  drawn  over  the  resolution  R (section  3.2).  Obviously,  Dec  and  Pee  are  both 
correlated  with  R : With  decreasing  R , the  decay  rate  and  the  phase  change  error  and  therefore 
numerical  dissipation  and  dispersion  increase. 

When  analysing  Dec  and  Pee  for  different  co,  it  is  not  astonishing,  that  cases  with  higher 
frequency  co  (smaller  k ) are  worse  resolved.  Nevertheless,  the  symbol  clouds  for  constant  co 
overlap  with  other  distributions,  because  the  incident  wave  angle  <p  has  a remarkable  effect  on 
the  resolution.  Values  of  R for  constant  co  can  even  differ  by  one  order  of  magnitude. 
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Using  logarithmic  scales  on  both  axes,  most  of  the  plotted  values  are  close  to  a straight  line 
( y = a ■ Rb  ).  The  parameters  a,b  can  be  obtained  by  the  least  squares  method: 

Declq  = 9.2  • R~25 
Pcelq  = 1.1  -R~2] 

From  linear  theory  for  a 2nd  order  CFD  code  on  a uniform  mesh  it  is  expected,  that  a 3rd  order 
relationship  exists  between  the  decay  rate  and  the  resolution  [6]  and  a 2nd  oder  relationship 
between  the  phase  change  error  and  the  resolution.  Better  agreement  with  the  3rd  order  decay  is 
achieved,  if  only  cases  B,C  and  D are  used  for  the  least  squares  fit  ( Declq  = 14.9  R~2  * ).  The 
results  of  cases  A might  be  influenced  by  the  convergence  criterion  (density  residual  Resp)  of 

the  code.  This  is  demonstrated  in  Figure  3.7,  where  simulations  with  good  convergence  criterion 
are  compared  with  simulations  with  worse  convergence  criterion:  For  higher  R , the  values  of 
Dec  and  Pee  deviate  significantly  from  the  least  squares  fit  of  the  good  convergence  criterion. 
For  a smaller  (better)  convergence  criterion  the  cases  with  lower  R are  almost  uninfluenced,  but 
the  values  of  Dec  and  Pee  for  higher  R come  closer  to  the  correlation  function. 

Dec  and  Pee  represent  the  accuracy  of  the  calculation  of  a wave  within  the  size  of  one  mesh 
cell  Ax . They  define  minimum  requirements  for  the  resolution  of  a pressure  wave.  As  an 
example:  if  the  decrease  of  the  pressure  amplitude  by  ldB  (0.1  dB)  due  to  numerical  dissipation 
after  200  mesh  cells  Ax  is  assumed  to  be  acceptable,  we  obtain  Dec  = 5 - 1 0~3  (5  • 1 0-4 ) and 
need  a resolution  of  the  wave  of  R - 20  (50). 

4.  Interaction  of  plane  waves  with  flat  plates 
4.1  Geometry  and  settings 


Figure  4.1:  Geometry  of  the  flat  plate  test  cases. 
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In  this  section,  the  transmission  and  reflection  of  plane  sound  waves  incident  upon  a single 
cascade  of  finite  plates  is  studied.  This  problem,  including  scattering,  has  been  solved  by  Koch 
[5]  by  means  of  the  finite  Wiener-Hopf  technique.  These  analytical  results  will  be  used  for 
comparison  with  computations  using  Lin3D. 

A single  cascade  of  parallel  flat  plates  of  length  l0  = lm , staggered  at  an  angle  fis  = 60°  and  at 
zero  angle  of  attack,  is  considered  and  is  shown  in  Figure  4.1.  (p  and  6 = <p~Ps  are  two 
different  definitions  of  the  angle  of  the  incident  wave.  For  the  computation,  the  regions  in  front 
of  the  cascade,  Lf,  and  behind  it,  Lh,  both  have  an  axial  extension  of  Lj  = Lh  = 0.5 m.  Two 

circular  frequencies  have  been  studied:  (0A  = 0). 25 Ttaj  as  a sub-resonant  example  and 

a>B  -7U. ij  as  a super-resonant  example,  where  scattered  modes  appear.  Only  cut-on  modes  are 

considered.  For  other  physical,  geometrical  and  numerical  parameters,  the  same  values  have 
been  chosen  as  for  the  wave  propagation  test  cases  in  section  3. 


Figure  4.2:  Unsteady  pressure  for  three  flat  plate  test  cases:  a)  coA , <p-  345° , b)  coA , 
<p  = 225° , c)  (oB , <p  = 225° . 


4.2  Sub-resonant  test  cases 

The  incident  waves  interact  with  the  blades,  when  passing  through  a blade  row.  If  the  blades  are 
perpendicular  to  the  wave  front,  the  wave  is  completely  transmitted.  For  blades  that  are  parallel 
to  wave  fronts,  the  waves  are  almost  completely  reflected  (Close  to  this  point,  full  reflection 
appears  for  a certain  angle  depending  on  the  Mach  number).  In  general,  the  incident  acoustic 
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wave  is  partially  reflected  and  partially  transmitted  through  the  blade  row.  For  the  sub-resonant 
test  case,  where  mode  scattering  does  not  appear,  two  snap-shots  of  the  pressure  field  are  shown 
in  Figures  4.2  a)  and  b). 


For  six  representative  cases,  some  parameters  can  be  found  in  Table  4.1.  Beside  the  incident 
wave  angles  0 , 6 , the  wave  numbers  k , kx,  ky  and  the  wave  lengths  X , Xx , Xy  are  listed: 


kv  = kx  X.an0 , X = . 


44 

'4+4 


. _ In  . _ 2k 


Waves  with  a cut-off  ratio  |£|>1  are  cut-on  (propagating  waves),  while  waves  are  cut-off 
(exponentially  decaying)  for  |£|  < 1 . 


Table  4.1:  Parameters  of  some  of  the  sub-resonant  flat  plate  test  cases. 
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Figure  4.3:  Variation  of  transmitted  (left)  and  reflected  pressures  (right)  with  incidence  angle 
for  the  sub-resonant  cases:  Lines:  analytical  solution  [5],  O Lin3D. 


The  cut-off  ratio  E,  is  also  shown  in  Table  4.1 : 
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Furthermore,  the  (analytic)  ratios  of  transmitted  pg  or  reflected  pg  wave  amplitudes  to  the 
amplitudes  of  the  incident  pressure  wave  p'0  are  listed  for  the  six  cases  in  Table  4.1. 


The  variations  of  these  ratios  | pg  j p'0  j and  |/>o incidence  angle  <p  are  drawn  in 
Figure  4.3.  A maximum  transmission  ratio  Ipjl  = 1 is  obtained  for  the  two  wave  directions, 


where  the  upstream  or  downstream  wave  fronts  are  perpendicular  to  the  blades  {<j>  = 60°,  240°). 

pd/|/?o|  = 0.  Analytical  and 


These  maxima  are  correlated  with  minima  of  the  reflection  ratio 
numerical  results  agree  well  with  each  other  within  the  whole  range  of  <p  ■ 


4.3  Super-resonant  test  cases 

Scattered  pressure  waves  ( m * 0 ) are  generated,  when  a pressure  wave  interacts  with  a cascade. 
For  some  of  the  cases  with  higher  circular  frequency  coB  = jmj  scattered  modes  can  be  cut-on: 

m = -1  is  cuton  for  1 0°  < <p  < 1 74° 

m = 1 is  cuton  for  223°  <<j>  < 294° 

Therefore  the  cases  with  this  frequency  are  called  super-resonant.  A snap-shot  of  the  pressure 
field  of  one  super-resonant  case  is  shown  in  Figure  4.2  c).  The  variation  of  the  ratios 

\prm  |/|/7/n  | with  incidence  angle  <fl  is  drawn  in  Figure  4.4.  Here  again,  analytical  and  numerical 
results  agree  well  with  each  other  within  the  whole  range  of  0 . 


5.  Summary 

The  application  of  the  time-linearized  Euler  method  Lin3D  to  2D  wave  propagation  problems  in 
homogeneous  flow  shows  the  dependency  of  numerical  errors  with  the  mesh  resolution.  For 
skewed  meshes,  the  minimum  number  of  mesh  diagonals  between  two  wave  fronts  is  an 
appropriate  measure  for  the  mesh  resolution.  Numerical  errors  o"f  the  discretization  scheme  in 
form  of  numerical,  dissipation  and  numerical  dispersion  increase  for  increasing  frequency  at 
constant  mesh.  For  constant  frequency,  the  resolution  of  a wave  depends  on  the  angle  of  the 
incident  wave  and  varies  within  one  order  of  magnitude.  The  angle  of  the  incident  wave  has 
therefore  a significant  effect  on  numerical  errors.  Correlations,  calibrated  with  the  simulations, 
can  be  used  further  to  check  the  mesh  size  before  performing  simulations  with  Lin3D. 

In  a second  study,  the  transmission  and  reflection  of  plane  sound  waves  incident  upon  a single 
cascade  of  finite  plates  has  been  calculated  with  Lin3D.  Computed  results  have  been  compared 
with  an  analytical  method  of  Koch  [5],  It  has  been  looked  at  sub-resonant  test  cases  and  at 
super-resonant  cases,  where  cut-on  scattered  modes  appear.  The  computed  ratios  of  transmitted 
or  reflected  to  incident  pressure  wave  amplitude  agree  well  with  the  analytical  solution.  Even 
the  ratios  for  scattered  modes  can  be  predicted  accurately  with  Lin3D. 

As  a next  step,  the  simulation  of  more  realistic  test  cases  with  finite  blade  thickness,  or  3D  flow 
is  envisaged. 


Figure  4.4:  Variation  of  transmitted  (top)  and  reflected  pressures  (bottom)  with  incidence 
angle  for  the  super-resonant  cases:  Lines:  analytical  solution  [5],  O Lin3D. 
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